import sys
from numpy import sqrt, zeros, cos, pi, log, arctan
import matplotlib.pyplot as plt
from scipy.linalg import solve
from scipy.integrate import quad
from scipy.special import yn

class BEMSolver:
    """
    Base class for constant boundary element method solver.
    
    
    xb - array of the x-coordinates of the boundary points
    yb - array of the y-coordinates of the boundary points
    
    """

    def __init__(self, xb, yb):
        self.elements = {}

        self.xb = xb
        self.yb = yb

        self.number_of_boundary_nodes = self.xb.size
        self.number_of_boundary_elements = self.number_of_boundary_nodes - 1

        self.L = zeros(self.number_of_boundary_elements)
        self.nx = zeros(self.number_of_boundary_elements)
        self.ny = zeros(self.number_of_boundary_elements)
        self.xm = zeros(self.number_of_boundary_elements)
        self.ym = zeros(self.number_of_boundary_elements)

        self.bctypes = zeros(self.number_of_boundary_elements)
        self.bcvalues = zeros(self.number_of_boundary_elements)

        for i in range(self.xb.size-1):
            self.xm[i] = 0.5 * (self.xb[i] + self.xb[i+1])
            self.ym[i] = 0.5 * (self.yb[i] + self.yb[i+1])
            self.L[i] = sqrt((self.xb[i+1] - self.xb[i]) ** 2.0 \
            + (self.yb[i+1] - self.yb[i]) ** 2.0)
            self.nx[i] = (self.yb[i+1] - self.yb[i]) / self.L[i]
            self.ny[i] = (self.xb[i] - self.xb[i+1]) / self.L[i]

        self.phi = zeros(self.number_of_boundary_nodes)
        self.dphi = zeros(self.number_of_boundary_nodes)

    def draw_boundary(self):
        """
        Create a matplotlib figure that illustrates the
        specified boundary points.
        """

        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.plot(self.xb, self.yb, 'o')
        plt.grid()
        plt.show()

    def cpf(self, xi, eta, x, y, nx, ny, L):
        raise NotImplementedError

    def solve(self):

        N = self.number_of_boundary_elements
        B = zeros(N)
        A = zeros((N, N), float)
        Z = zeros(N)

        for m in range(N):
            B[m] = 0.0

            for k in range(N):
                PF1, PF2 = self.cpf(self.xm[m], self.ym[m], self.xb[k],
                                    self.yb[k], self.nx[k], self.ny[k],
                                    self.L[k])
                F1 = PF1 / pi
                F2 = PF2 / pi

                if self.bctypes[k] == 0:
                    A[m, k] = -F1
                    B[m] += self.bcvalues[k] * (-F2 + 0.5 * delta(k, m))
                else:
                    A[m, k] = F2 - 0.5 * delta(m, k)
                    B[m] += self.bcvalues[k] * F1

        Z = solve(A, B)

        for m in range(N):
            if self.bctypes[m] == 0:
                self.phi[m] = self.bcvalues[m]
                self.dphi[m] = Z[m]
            else:
                self.phi[m] = Z[m]
                self.dphi[m] = self.bcvalues[m]

    def evaluate_phi_at_interior_point(self, xi, eta):
        sum_ = 0.0

        for i in range(self.number_of_boundary_elements):
            xb = self.xb[i]
            yb = self.yb[i]
            nx = self.nx[i]
            ny = self.ny[i]
            L = self.L[i]
            PF1, PF2 = self.cpf(xi, eta, xb, yb, nx, ny, L)
            sum_ += self.phi[i] * PF2 - self.dphi[i] * PF1

        return sum_ / pi


def delta(m, n):
    """
    Kroneckers delta
    """
    if m == n:
        return 1.0
    else:
        return 0.0

class TwoDimensionalLaplaceBEMSolver(BEMSolver):

    def cpf(self, xi, eta, x, y, nx, ny, L):
        A = L ** 2.0
        B = 2.0 * L * (-ny * (x - xi) + nx * (y - eta))
        E = (x - xi) ** 2 + (y - eta) ** 2
        D = sqrt(abs(4.0 * A * E - B ** 2))

        if D < 1e-12:
            PF1 = 0.5 * L * (log(L) \
            + (1.0 + 0.5 * B / A) * log(abs(1.0 + 0.5 * B / A)) \
            - 0.5 * B / A * log(abs(0.5 * B / A)) - 1.0)
            PF2 = 0.0

        else:
            PF1 = 0.25 * L * (2.0 * (log(L) - 1.0) \
            - 0.5 * B / A * log(abs(E / A)) \
            + (1.0 + 0.5 * B / A) * log(abs(1.0 + B / A + E / A))
            + D / A * (arctan((2.0 * A + B) / D) - arctan(B / D)))

            PF2 = L * (nx * (x - xi) + ny * (y - eta)) / D \
            * (arctan((2.0 * A + B) / D) - arctan(B / D))

        return PF1, PF2


